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Many-electron systems confined at substantial finite temperatures and densities present a ma- 
jor challenge to density functional theory. Comparatively little is known about the free-energy 
behavior of such systems for temperatures and pressures of interest, for example, in the study of 
warm dense matter. As a result, development of approximate free-energy functionals is faced with 
difficulties of assessment and calibration. Here we address, in part, this need for detailed results 
on well-characterized systems. We present results on a comparatively simple, well-defined, but 
computationally feasible model, namely the thermal Hartree-Fock approximation applied to eight 
one-electron atoms with nuclei at fixed, arbitrary positions in a hard-walled box. We discuss the 
main technical tasks (defining a suitable basis and evaluation of the required matrix elements) and 
discuss the physics which emerges from the calculations. 



I. INTRODUCTION 

Warm dense matter (WDM) is encountered in systems 
as diverse as the interiors of giant planets^ and in the 
pathway to inertial confinement fusion 3,4 . WDM is chal- 
lenging to theory and simulation because it occurs in- 
conveniently, for theory, between the comparatively well- 
studied plasma and condensed matter regimes. Both the 
Coulomb coupling parameter T := Q 2 /(t- s £;bT) and elec- 
tron degeneracy parameter O := fe^T/e^ are approx- 
imately unity for WDM. (Q = relevant charge, r s = 
Wigner radius, Bp = electron Fermi energy, T = temper- 
ature, fcs = Boltzmann constant.) A non-perturbative 
treatment therefore is required. 

Contemporary computations on WDM^— are domi- 
nated by use of the Kohn-Sham (KS) realization of ther- 
mal density functional theory (DFT)^~— to generate a 
potential surface for ionic motion (treated classically). 
The great majority of these calculations use approximate 
ground-state exchange-correlation (XC) functionals, E xc , 
with the temperature dependence of the XC free energy 
picked up implicitly from the T-dependence of the den- 
sity n(r,T). Though fruitful, this approach is not with- 
out potential difficulties, as is illustrated in Fig. 3 of Ref. 
0. Three issues are germane here. 

First, virtually nothing systematic is known about the 
implicit T-dependence of ground-state approximate E xc 
functionals (especially beyond the local density approxi- 
mation, LDA), whether they be constraint-based or em- 
pirical. Compared to the ground-state situation, there 
is only a small literature on explicitly T-dependent func- 
tionals, that is, XC free energy functionals, and essen- 
tially all of those studies are at the level of the LDA 2 ^— . 

Second, there is the computational burden of solving 
for the Kohn-Sham orbitals and eigenvalues. Since the 
computational load from the eigenvalue problem scales, 
in general, as order A r ( 3 rbitals , the growth in the number of 
non-negligibly occupied KS orbitals with increasing tem- 
perature is a clear computational bottleneck. See the re- 



marks, for example, in Sec. 4 of Ref. [l8|. For complicated 
systems, the same bottleneck is encountered in ground- 
state simulations which use the DFT Born-Oppcnhcimcr 
energy surface to drive the ionic dynamics. One result 
has been the emergence of active research on orbital-free 
DFT (OFDFT), that is, approximate functionals for the 
ingredients of the KS free energy, namely the KS kinetic 
energy (KE) T s , entropy S s , and XC free energy T xc or 
their ground-state counterparts. Almost all of this ef- 
fort has been for ground-state OFKE functionals 4 ^—. 
(Note that most of the OFKE literature invokes the KS 
separation of the KE in order to use existing E xc approx- 
imations consistently.) 

Third, the finite-temperature OFDFT work is domi- 
nated by variants on Thomas-Fermi-von Weizsacker the- 
ory; see for example Ref. 41 and references therein. That 
type of theory, however, is known (on both fundamental 
and computational grounds) to be no more than qualita- 
tively accurate in many circumstances relevant to WDM 
(e.g., chemical binding). Compared to the data-rich 
context for development of zero-temperature functionals, 
there is little to guide development and assessment of 
finite-T functionals. Similarly, compared to the T=0 sit- 
uation (or the very high T situation), not much is known 
about the accuracy of approximate finite-T OFDFT func- 
tionals beyond TFvW. An aim of the present study is to 
provide such data for the combination of a well-defined 
physical system with a well-defined approximation and 
its implementation. 



II. SYSTEM AND METHODOLOGY 

We shall treat a neutral system of Ni on atoms with 
nuclei fixed at arbitrary positions in a hard-walled three- 
dimensional rectangular box. The confined system allows 
systematic treatment of pressure effects at stipulated fi- 
nite temperature, hence is a small, treatable sample of 
WDM. For specificity in the initial work, we chose one- 
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electron atoms and Ni on < 8 for simplicity. Note the 
corresponding use of a small number of atoms in Ref. 
IT3 , Also note the considerable literature on spherically 
confined systems at T = K^Zr— . We have not found 
any work on lower-symmetry confinement of multi-atom 
systems at non-zero T. At non-zero T there is a large 
literature on average-atom methods, for example Refs. 
[18II53I and many others, but such methods do not provide 
the fiduciary data needed for functional development. 

As a many-electron problem, a clearly defined approx- 
imation is required. We choose the finite temperature 
Hartree-Fock (FTHF) scheme^^^, with issues of elec- 
tron correlation to be addressed in the future. FTHF 
provides a clear definition for the correlation free en- 
ergy which is consistent with the definition used in zero- 
temperature theory with approximate wave functions. 
(Of course, the definition of correlation energy in zero- 
temperature DFT differs, as it does in finite-temperature 
DFT, but that is no more a barrier here than there.) 
In addition to being well-defined in the grand ensemble, 
FTHF provides the advantage that its T = K limit is 
the lingua franca of molecular electronic structure inter- 
pretation. Use of FTHF therefore also provides at least a 
semi-quantitative framework for understanding chemical 
processes in WDM. 

The FTHF approximation is defined in the grand 
canonical ensemble by restricting the relevant traces to 
states which are single Slater determinant o 23 i 54 ' 55 . The 
result is an upper bound to the free energy Ffthf > J~. 
Standard thermodynamic relationships for the grand en- 
semble follow. The FTHF Euler equations to be solved 
(in unrestricted form) are^ 

em{r) = (y-\^ 2 + U ion (r) \ ^(r) 
j 

with Ui on the ion-electron interaction potential and Oi 
the spin label; the sums are over all spin orbitals. Un- 
less indicated otherwise we use hartree atomic units 
(h = m e = e = 1; energy is then in hartrees, 1 hartree = 
27.2116 eV, and lengths are in bohrs, 1 bohr = 0.52918 
angstrom). The spin orbitals (eigenstates, (pi) have 
Fermi-Dirac thermal occupations 

fi= (l + e^-^y 1 , N = J2fi, (2) 

i 

where /3 = l/fc^T and [i is the electron chemical poten- 
tial. For a specified value of N, which is a grand ensemble 
average, /i must be determined. These equations, along 
with specification of the nuclear sites and imposition of 
hard-wall boundary conditions, completely describe the 
problem. 



The FTHF free energy and entropy are given by 

Ffthf = ^ fi£% - g X! h (^*>j _ ~ TSfthf , 

(3) 

Sfthf = -k B fi M/0 + (1 - fi) Ml - fi) > ( 4 ) 

i 

with conventional definitions of J and K: 

(e) 

Equations ([3]) and (j4]) clearly reduce to conventional 
ground-state Hartree-Fock expressions in the zero- 
temperature limit. 

A. Gaussian Basis Sets 

Solution of self-consistent field equations such as (fi]) 
via Gaussian-type-orbital (GTO) basis methods is the 
standard procedure in modern computational codes for 
molecular systems. First introduced by Boys^— , such 
basis sets automatically satisfy the free-molecule bound- 
ary condition that the orbitals vanish at infinity. For 
a hard-wall confined system, the basis functions must 
vanish at the boundary, so standard molecular GTO 
matrix-element expressions are inapplicable. This fact 
illustrates the most critical technical issue for implemen- 
tation, namely, to find a basis that satisfies the bound- 
ary conditions yet allows for an efficient enough evalua- 
tion of the two-electron integrals to be computationally 
tractable on reasonable resources. A second, closely re- 
lated technical issue is that the high temperature also dic- 
tates what is "efficient enough" , in that the basis must be 
large enough and flexible enough to represent a sufficient 
number of thermally occupied higher-energy orbitals of 
the system, hence to represent the density and free en- 
ergy accurately. 

Those considerations eliminate several seemingly plau- 
sible options for a basis. For example, a real-space finite- 
difference/element scheme, while suitable for a DFT 
calculation or a Hartree-Fock calculation on a free di- 
atomic molecule (for which curvilinear coordinates can 
be exploited^), is far too expensive for the present case 
because of the number of matrix elements to be calcu- 
lated. Another example is sine functions. They also sat- 
isfy the hard-wall boundary condition, but an adequate 
description of the rapidly varying electronic distribution 
near the nuclei requires prohibitively many matrix el- 
ements in our multi-center problems. So we chose an 
adaptation of standard GTO methods which uses mod- 
ified Gaussians that meet the boundary condition, yet 
retain enough efficiency to complete the calculation. 
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The various ways to force a GTO to zero at the bound- 
ing planes of a rectangular box can have great impact 
upon the efficiency of the matrix element calculations. 
Compared with familiar practice for free molecules, in 
general the confined case requires more primitive GTOs 
for each contracted one. More importantly, the finite 
integration volume makes it impossible to achieve com- 
pletely analytic calculation of the two-electron integrals, 
which is, of course, precisely the category in which com- 
putational efficiency is most needed. We have addressed 
this issue by using truncated Gaussians as described in 
the next section. 



B. Truncated Gaussians 

The rectangular box makes Cartesian GTOs a conve- 
nient choice, because each primitive function then is sep- 
arable into Cartesian factors which are simple ID func- 
tions. Consider the Cartesian factor 

g n (x) = (x- x c ) n e- a{x ' x ^ 2 (7) 

To force this function to zero at the box boundaries x = 
and x = L x , we subtract a constant equal to the function 
value at each end. When x c is not at the box center, the 
value to be subtracted differs for the two ends, so we split 
the function into two pieces, make the two subtractions, 
and scale the two pieces such that the resulting func- 
tion is continuous. Each unnormalized Cartesian factor 
becomes 

gl ox (x) = a (g n (x) - A ) < x < x c 

= a L (g n {x)-A L J x c <x<L x (8) 

with A = 5™(0), A Lx = g n (L x ). We call this the trun- 
cated GTO (tGTO) basis. 

Two technical issues remain. The tGTO functions may 
not have continuous derivatives, so proper evaluation of 
the kinetic energy matrix elements requires attention. 
Appendix A shows that nothing untoward happens and 
that the kinetic energy is simply a sum of piecewise con- 
tributions, except for p-type functions which have a sim- 
ple correction term. Second is the matter of evaluating 
two-electron matrix elements. In Appendix B, we show 
that this task reduces to computing finite-range integrals 
of products of Gaussians and error functions. At this 
juncture, we are doing those via Gauss-Legendre quadra- 
ture. We also note that the tGTO basis is simpler than 
the smoothly cut-off floating spherical GTO basis of Ref. 
[5ll in the sense that their cut-off introduces a mixing of 
two symmetry types (e.g. s, d) in each basis function. 

Note that, so far, we have implemented only 
the restricted Hartree-Fock approximation (RHF; non- 
spin-polarized in DFT language; closed-shell or spin- 
compensated in quantum chemistry language). 



III. RESULTS 
A. Zero Temperature 

Two simple ground-state test cases, the H atom and 
the H2 molecule, illustrate the system behavior with in- 
creasing confinement (decreasing volume) as well as the 
correctness of our implementation. The confined-system 
energies should be above the ground-state energies (in 
the basis selected) of the corresponding free systems and 
approach those free-system energies (and bond length for 
the molecule) in the limit of large box volume. 

For the hydrogen atom in the center of a cubical box, 
Fig.[T](a) shows the ground-state energy for confinements 
L < 10 bohr. (The basis exponents shown in Table Q] 
their selection is described in the next section.) The cal- 
culations were carried out to L — 30 bohr. Beyond 10 
bohr there was negligible difference in the energy with 
respect to the free-atom energy in the same basis, pre- 
cisely as expected. At L = 30 bohr the ground-state 
energy is identical with the free-atom GTO calculation 
using the same basis, namely —0.498476 hartree, vali- 
dating our overlap, kinetic, and nuclear energy integral 
calculations. A fit to the energies as a function of cube 
edge L 

E(L)=a/L 2 + b/L + c , < L < 10 bohr (9) 

yields a — 14.5733 hartree bohr 2 , b = —3.82369, hartree 
bohr, c = —0.238258 hartree. Though this fit does not 
have the correct infinite-size limit, it is the best fit of 
this simple form for < L < 10 bohr. ^From this fit the 
pressure, p = —dE/dV (V — L 3 ) can be calculated; see 
Fig- [11(b). 

The effects of spherical versus cubic confinement can 
be assessed easily by comparison of the two confinement 
types both for bounding volumes and equal volumes. For 
bounding volumes, the circumscribed-sphere R = y / 3L/2 
and inscribed-sphere R = L/2 results from Ref. |59j (also 
see Ref. HI) can be compared to our cube results. Figure 
[5] shows that the energy of the cubically confined system 
is bounded by that of the two spherical systems, as ex- 
pected. Shape effects of confinement are shown in Fig. 
[3J Though the spherical system has a lower energy than 
that of the equal volume cubical system, they are quite 
close until the cube is smaller than L s» 4 bohr, where 
a significant indication of shape dependence begins. At 
L = 1 bohr the difference in energy is over 25%, as shown 
in Fig. 131 (b). A simpler basis of six s-type tGTO with 
exponents [0.15,0.3,0.6,1.2,2.4,4.8 bohr -2 ] reproduces the 
ground-state energies for L > 1.5 bohr. 

Next consider the confined H2 molecule at zero tem- 
perature. Here we used the simpler six-tGTO basis just 
given. For a large cube, the new tGTO confined-box 
computations again should conform to known results for 
the integrals and produce essentially the energy vs. bond 
length curve for the free molecule. With L — 30 bohr 
and the molecule centered in the cube and aligned along 
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FIG. 1: (a) Ground-state energy of an H atom in the center 
of a cube with edge length L. Fit is for function f(L) — 
a/L 2 + b/L + c. (b) Pressure calculated from the energy fit. 



FIG. 3: (a) Comparison of energy of an H atom in a cube 
with that in a sphere of equal volume, (b) The difference of 
the two ground-state energies. 




FIG. 2: Comparison of ground-state energy of an H atom in 
a cube with its energy within two bounding spheres. 



the body diagonal, we get energies shown as points in 
Fig. |U These agree completely with the values of the 
free GTO calculation, which are shown in that figure as 
the "Free" curve. The minimum is at 1.383 bohr, sat- 
isfactorily close to the free-molecule RHF value of 1.385 
bohr from a 6-31G** basis calculation—. Conversely for 
fixed R (at 1.4 bohr) and decreasing L (to 4 bohr), the 
ground-state energy behavior is shown in Fig.O The on- 
set of confinement effects becomes visible in the vicinity 
of L s» 11 bohr. Optimization of the bond length R at 
L = 5 bohr is shown in Fig. [5] The total energy is, of 




-1.125 1 1 1 1 1 1 1 1 1 

1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 

Separation (bohr) 

FIG. 4: Energy versus atomic separation for the H2 molecule 
centered in a 30 bohr cube. The ground-state agrees exactly 
with the free GTO calculation using the same primitives as a 
basis. 



course, higher than for the larger box, and the optimal 
R is shifted down to 1.178 bohr from 1.383 bohr. 

Following this method, we obtain the optimized R and 
energy for decreasing L. A function fit to the energy 
similar to that used earlier yields R as a function of the 
pressure, as shown in Fig. [7] 
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L (bohr) 




FIG. 5: Ground-state energy for H2 with bond length 1.4 
bohr confined in cube of edge length L. 



FIG. 8: Lowest five eigenstate energies of H as a function of 
cube edge L. Energy values are from Table U 
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FIG. 6: Energy vs. bond length for the H2 molecule centered 
in cubes L — 30 and L = 5 bohr. The minima are at 1.383 
and 1.178 bohr respectively. 



B. Finite Temperature 

For finite-temperature calculations, at least a single p 
orbital needs to be included, and as many orbitals as 
feasible should be available to represent the fractionally 
occupied levels which become increasingly important as 
T is increased. For all the following calculations, except 
where noted, the basis consists of seven s-type primi- 
tive GTOs and one p x , one p y , and one p z GTO, with 
the p-GTO exponents equal. An elementary exponent 




100 200 300 400 500 600 
Pressure (Gpa) 

FIG. 7: Optimized bond length, R, for the H2 molecule vs. 
pressure. 



optimization was done as follows. A set of five s-type 
exponents was picked and held fixed: [0.1,0.2,0.4,0.8,1.6 
bohr -2 ]. The sixth exponent was optimized to to mini- 
mize E\ s for the single atom centered in a cube of specified 
L. With those six exponents fixed, the seventh s-GTO 
exponent was used to minimize e 2 s- With those seven 
fixed, the p exponent was used to minimize S2 V - Addi- 
tionally, for small L the smallest exponents produce or- 
bitals that are so similarly flat that an approximate linear 
dependence exists and diagonalization fails. When this 
happens the smallest exponent is replaced by extending 
the even-tempered exponent series to larger values. For 
example at L = 4 the 0.1 exponent is replaced with 3.2, 
at L — 2 the 0.2 exponent is also replaced with 6.4. 

This procedure keeps the ratio of the effective length 
(1/y/ci) of the most diffuse function to the edge length L 
at 0.75-0.79 for L < 3, with smaller ratios for larger L. 

The optimization was done for each L. Table U shows 
the resulting exponents and energies for the orbitals that 
were optimized. 

With this volume-dependent optimized basis, calcula- 
tions were done for a single atom at the center of a cube 
with 1 < L < 15 bohr. The orbital energies of the five 
lowest states (Is, 2p Xl 2p yi 2p z , 2s) are plotted in Fig. [8] 
Notice the inversion of ordering (2p below 2s) that is a re- 
sult primarily of the confinement. To address finite tem- 
perature for this single-electron system, the one-electron 
levels were populated according to the Fermi-Dirac dis- 
tribution. Observe that the one-electron Hamiltonian is 
independent of density, so the one-electron orbitals and 
eigen-energies are independent of occupancy, even though 
the density and total electronic energy are not. The left- 
hand panel in Fig. |9] shows the resulting total energy as 
a function of L for four values of T, while the right-hand 
panel shows the free energy. The weak minimum in total 
energy in the vicinity of 6 bohr at T=50 kK appears to be 
a confinement effect. We have found a similar minimum 
at about the same volume by doing a Fermi-Dirac popu- 
lation of the high-precision eigenvalues of the spherically 
confined H atom given in Ref. [59T 

FigurefTUlshows the contributions to the free energy for 
the single atom as a function of T for four cube volumes 
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FIG. 9: (a) Single H atom total energy, E = Tfthf + TSfthf, and (b) free energy, Tfthf as a function of cube edge L, 
with the one-electron levels populated according to the Fermi distribution. 
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FIG. 10: Energy components of a cubically confined single H atom vs. temperature, (a) kinetic energy, (b) electron-nuclear 
energy, (c) total energy, (d) free energy, (e) entropic energy [(d)-(c)]. Key for all is shown in (e). 
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L 


fixed 


Is 


2s 


2p 


Els 


£■23 


£2p 


1.0 


1.6, 3.2, 6.4, 12.8, 25.6 


179.1 


244 (NA) 


4.01 


10.518 


50.5898 


26.683 


1.5 


0.8, 1.6, 3.2, 6.4, 12.8 


92.15 


244 (NA) 


1.88 


3.6738 


21.3545 


11.154 


2 


0.4, 0.8, 1.6, 3.4, 6.8 


48 


250 


1.115 


1.48471 


11.3649 


5.8753 


3 


0.2, 0.4, 0.8, 1.6, 3.4 


24.5 


175 


0.545 


0.11385 


4.47073 


2.25356 


4 


0.2, 0.4, 0.8, 1.6, 3.4 


21.2 


7.0 


0.34 


-0.268848 


2.18313 


1.06314 


5 


0.1, 0.2, 0.4, 0.8, 1.6 


11.2 


3.7 


0.235 


-0.40474 


1.18372 


0.547955 


6 


11 


10.4 


2.5 


0.18 


-0.458898 


0.675591 


0.288794 


7 


11 


17.7 


2.25 


0.15 


-0.481704 


0.389716 


0.145105 


8 


» 


10.2 


0.195 


0.122 


-0.491112 


0.217062 


0.0591604 


9 




10.2 


0.07 


0.103 


-0.495291 


0.10651 


0.00450443 


10 


» 


10.2 


0.0365 


0.088 


-0.497104 


0.0327616 


-0.0321549 


11 




10.1 


0.023 


0.076 


—0.497889 


—0.0179536 


—0.0578267 


12 


11 


10.1 


0.018 


0.0665 


-0.498227 


-0.0535872 


-0.0764093 


13 


11 


10.1 


0.015 


0.059 


-0.498372 


-0.0789901 


-0.0901715 


14 


11 


10.05 


0.014 


0.053 


-0.498435 


-0.0972631 


-0.100500 


15 


11 


10.1 


0.014 


0.048 


-0.498461 


-0.108291 


-0.108291 



TABLE I: Optimized exponents and orbital energies for an H atom at the center of a cube of edge length L. 2p refers to the 
triply degenerate p x , p y , and p z states. For L < 2, the 2s level could not be optimized beyond the first 6 exponents. 



(L = 2, 3, 7, 15 bohr). At L = 7 bohr, the KE is flat with 
T at almost its T = K value. The T = 50 kK nuclear- 
electron attraction E^ ei however, is much stronger for 
L = 7 bohr than for L = 15 bohr. By L = 3 bohr, the 
KE and Ej^ e are roughly equal in magnitude. 

Figure [TOl also shows that the KE for the L = 15 bohr 
system falls with increasing temperature. Though this 
might seem odd, it is as it should be from virial theorem 
arguments for the free atom. The 2s KE is one-fourth 
the Is valued. Finite temperature population of the 2s 
and depopulation of the Is therefore reduces the KE with 
respect to its T = K value. 

Next we turn to the system of eight H atoms. We 
examined a symmetric configuration in which the eight 
atoms were situated at the corners of a smaller cube, 
edge-length L/2, centered within the hard- wall cube, 
edge-length L. The basis used was ten s-type GTOs cen- 
tered on each atom. Strict s-type symmetry is broken, of 
course, by enforcement of the hard-wall boundary condi- 
tions. An even-tempered set of exponents also was used 
in this case: [0.2,0.4,0.8,1.6,3.2,6.4,12.8,25.2,50.4,100.8 
bohr -2 ]. As a test, the calculations were redone with two 
fewer basis functions per atom; the exponents [50.4,100.8 
bohr -2 ] were removed. The two calculations agreed to 2 
millihartree in total and component energies up to 200 
kK. Matrix elements are calculated only once and stored, 
after which fully self-consistent calculations may be done 
at many temperatures. 

Figure QT] shows the total energy E = Ffthf + 
TSfthf as a function of L for various temperatures, as 
well as the free energy Ffthf itself. The nuclear-nuclear 
repulsion energy is included (constant with respect to 
temperature, it varies with L). Figure [T2l shows various 
components of the free energy (electron-nuclear Coulomb 
energy, electron-electron Coulomb energy including ex- 
change, kinetic energy, and entropic energy). Also shown 
are the total energy and free energy, with the difference of 



these two being the entropic energy. Again, the nuclear- 
nuclear repulsion energy is included in these two plots. 
It has values of 9.118, 7.598, 5.699, 4.559 hartree for L = 
5, 6, 8, 10 bohr, respectively. The energies are shown 
as a function of temperature < T < 250kK for four 
cube sizes, L — 5, 6, 8, and 10 bohr. The compara- 
tively flat plateau up to roughly T = 15 kK is a direct 
consequence of Fermi-Dirac level filling. Up to about 25 
kK, the interval between the highest occupied molecular 
orbital (HOMO) at zero temperature and the lowest un- 
occupied molecular orbital (LUMO) is roughly constant 
at 0.5 hartree. Therefore the filling ratio of those two is 
roughly exp(— 13/1.3) ~ 5 x 10~ 5 or smaller up to about 
T = 15kK. We return to this matter below. 



C. Comparison with approximate functionals 

Orbital-free treatment of WDM has been dominated, 
not surprisingly, by local density approximations for both 
the KE and exchange contributions to the free energy. 
For the KE, the choice is physically motivated by the fact 
that the high- pressure and/or high-temperature limit for 
a WDM system is Thomas-Fermi. This leads to finite- 
temperature Thomas-Fermi (FTTF^l % = J t dr by 
itself, or with some fraction of von Weizsacker contribu- 
tion (in its zero-temperature form) Tw — J T w dr, with 

^ = ^ph/tm (io) 

n = r2fl372 7 V2 W (11) 



where the I are Fermi integrals. Note a parametrized 
form— may be used to eliminate fi between tq and n. 
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FIG. 11: Zero and finite-temperature total energy and free energy for eight H atoms in different sized cubes: (a) E 
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FIG. 12: Self-consistent FTHF free energy contributions for eight confined H atoms in a cubic array as a function of T for four 
different cube edge lengths L. Key for all is given in (d). 
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FIG. 13: FTHF kinetic energy compared with finite- 
temperature Thomas-Fermi KE and two forms of von 
Weizsacker augmentation of FTTF for L — 6 bohr. See text 
for details. 



In Fig. [13] we compare the FTHF KE as a function of 
T with the FTTF KE alone and with it supplemented by 
both the full Tw and (l/9)Tw for L = 6 bohr. The latter 
three functionals were evaluated with the FTHF density, 
hence are non-self-consistent. As is known at T = 0, pure 
Thomas-Fermi underestimates the KE while addition of 
the full Tw overestimates it. None of the three is close to 
quantitative agreement with Tfthf- Moreover, FTTF 
and FTTF augmented with (1/9)7W can be ruled out 
from the T = 0K comparison, since the exact (fully cor- 
related) KE must be above Tfthf (from a virial theorem 
argument). 

For the exchange contribution to the free energy, the 
use of ground-state functionals is common (recall In- 
troduction), with the LDA being dominant. Fig. Q3] 
shows the FTHF exchange contribution to the free en- 
ergy J~ x, fthf in comparison with the exchange free- 
energy generated by ground-state LDA functional and 
with the Perrot and Dharma-wardana parametrization 
for the temperature-dependent LDA functional- 8 . Again, 
this is for L = 6 bohr and the other functionals are evalu- 
ated with the FTHF density. (Note that such "post-scf" 
evaluation is fairly common in assessment of newly devel- 
oped ground state functionals. See, for example, Ref.l68h 
Here one sees a marked difference: the ground-state func- 
tional fails completely while the temperature-dependent 
functional has at least semi-quantitatively correct tem- 
perature dependence. 

We may also make some semi-quantitative comparison 
with a more widely used model for extended systems at 
substantial T. In Fig.[TS]we show the internal energy per 
atom of the eight hydrogen atoms in the cubic symme- 
try arrangement at L — 7 bohr (corresponds to average 
r s = 2.17), with that of the average atom DFT calcula- 
tion of Dharma-wardana and Perrot— at r s = 2. Their 
system includes DFT exchange and correlation whereas 
we have pure Hartree-Fock exchange. Additional energy 
differnces are due to the different boundary conditions. 
Decomposing the near parallel temperature dependence 
into those components, however, is a task outside the 



FIG. 14: FTHF exchange free energy ground-state LDA and 
Perrot and Dharma-wardana temperature-dependent LDA 
(L = 6 bohr). See text for details. 
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FIG. 15: FTHF internal energy per atom for L = 7 bohr, as 
compared with r s = 2 hydrogen plasma average atom DFT 
calculation. 

scope of the present work. We can already see one as- 
pect. Though the kinetic energy is the major contributor 
to the change in total energy as a function of tempera- 
ture (recall Fig. I12p. the change due to electron-electron 
interaction, including exchange, is at least a third that 
of the kinetic energy. 



IV. DISCUSSION AND CONCLUSIONS 

Comparison, evaluation, and betterment of function- 
als for WDM simulations is the long-term motive of this 
work. As just shown, even at this initial state of develop- 
ment (Hartree-Fock, small particle number, no molecu- 
lar dynamics), the approach gives insight regarding that 
task. There are some specific issues worth discussion also. 



A. One-particle spectrum effects 

It is well known that zero-temperature HF calcula- 
tions over-estimate both band gaps in solids and the 
so-called HOMO-LUMO (highest occupied and lowest 
empty molecular orbital respectively) gaps of finite sys- 
tems. This happens because the occupied HF orbitals 
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are free of self-interaction, but the unoccupied ones are 
not. As temperature increases in the FTHF scheme, how- 
ever, levels unoccupied at T = K become increasingly 
occupied and shifted. In the self-consistent solution of 
the HF problem in a basis, those levels then contribute 
to the Hamiltonian matrix ("Fock matrix" in quantum 
chemistry terminology). Thus there are two questions 
to address. What is the extent to which FTHF ex- 
hibits HOMO-LUMO gap over-estimation similar to that 
of ground-state HF? At what fractional occupancy is an 
energy level changed from an overestimated virtual to a 
more properly estimated partially occupied one? 

For specificity, we consider the eigenspectrum of a sin- 
gle, moderately compressed eight-atom cube with L = 6 
bohr. Figure [16] shows the Fermi distribution of the 
single-particle energies for four temperatures. Note that 
these distributions are for one spin. Also keep in mind 
that the cubical symmetry causes the lowest four one- 
electron orbitals to group as singly degenerate and triply 
degenerate (alg, tlu in crystal field notation). At zero 
temperature therefore, only two points are shown with 
unit occupancy, but the higher energy point corresponds 
to three degenerate HOMO states (indexed as 2, 3, 4). 
The LUMO is the degenerate states 5, 6, and 7, with the 
singly degenerate state 8 above them (again as would 
be expected from a cubical crystal field). For simplic- 
ity of discussion, the HOMO and LUMO (at T = 0) 
are labeled £4 and £5 respectively. One can see that the 
spacing between £4 and £5 decreases with temperature. 
This difference is shown directly in Fig. [18] along with 
the occupation number for the £5 level. The continuous 
curves in Fig. ll6l show the Fermi function with calculated 
chemical potential fi. The discrete points mark the input 
energies to the calculation of /.i. Then in Fig. [TTl /i is 
plotted as a function of T. Observe that the chemical 
potential is nearly mid-way between £4 and £5 up to just 
below 50 kK. This behavior is exact at zero temperature. 

Those total energy and eigenspectrum results together 
resolve the matter of the behavior of what would be vir- 
tual states at zero temperature in FTHF. First, exami- 
nation of the kinetic and total energy plots for all box 
sizes shows that there is a change in the form of the tem- 
perature dependence at roughly 20 kK. That change is 
complemented by the change in the spacing between the 
£4 and £5 levels. They are essentially static for lower 
temperatures, then begin to change abruptly well below 
50 kK, and then change more moderately at higher tem- 
peratures. Thus, above 50 kK, states corresponding to 
zero-temperature virtuals are sufficiently incorporated in 
the interaction terms to make a material modification 
of T = behavior. However, as the temperature is de- 
creased below roughly 50 kK, the FTHF Coulomb and 
exchange terms increasingly are dominated by the T = 
occupied levels, which therefore keeps the lightly occu- 
pied higher energy levels artificially high. This is not 
a basis issue, but an issue with discrete eigenstates. In 
a solid such as jellium or a metal with a continuum of 
states, this should not be an issue, but for a system with 
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energy gaps the issue remains. 



B. Other findings and considerations 

The preceding discussion about depopulation and re- 
population of levels relative to the ground-state HF 
HOMO-LUMO gap illustrates a broader challenge for 
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FIG. 18: Energy difference between the fourth and fifth en- 
ergy levels along with the occupation number of the fifth en- 
ergy level. 
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the construction of approximate functionals for the vari- 
ous contributions to the free energy. Whether explicit or 
implicit, such approximate functionals correspond to re- 
stricting the required traces to specific classes (or sub-sets 
of classes) of state functions. The consequence of such a 
restriction is to incorportate the spectral properties of 
that class into the approximate functionals. For exam- 
ple, in FTHF that class is single Slater determinants con- 
structed with respect to ground-state HF minimization. 
Although we have not attempted its construction here, 
in principle there is an FTHF free-energy functional. It 
would have exactly the same problem with a plateau in 
its T-dependence as we have found here. 

The small number of particles is another issue. For 
sufficiently large numbers, all standard ensembles (grand 
canonical, canonical, micro-canonical) give the same 
thermodynamics. Fluctuations characteristic of small 
particle counts can degrade that relationship. Such issues 
are endemic to any computational study, especially when 
a large temperature and pressure domain such as charac- 
terizes WDM is involved. In that sense, the present study 
is not much worse off than most modern WDM simula- 
tions with small numbers of electrons^. At the least, 
we have an even-handed comparison of different meth- 
ods (e.g., the comparison of functionals given above) for 
a given number of electrons and of ions. The main issue 
regarding particle count is computational cost. 

Clearly we have shown that the tGTO basis is feasi- 
ble and effective. As is typical of GTO basis methods, 
the computational cost is essentially entirely in the cpu 
time for the calculation of the two-electron integrals. The 
eight- atom systems described are calculated with 64 or 
80 total basis functions, making diagonalization trivial. 
A simple double array of all A^ 4 two-electron integrals 
only occupies 128 or 312.5 MB respectively, storable in 
memory. In practice we calculated N 4 /X integrals, with 
A =7.60, 7.68 and not 8 due to the looping procedure we 
used. While for the cubic system discussed here, this cal- 
culation requirement could have been reduced further by 
exploiting symmetry, we need the capability to explore 
other, lower symmetry geometries. Note however that 
if the box size or the atomic positions are changed, all 
integrals affected by the change must be recalculated. 

On a modern desktop processor (Intel Core i5 650 at 
3.2 Ghz) one two-electron integral can be calculated in 
about 29.1 ms. So the times to calculate all integrals 
for a 64- or 80-orbital basis would be 17.84 hr and 43.13 
hr. The integrals are calculated independently, so can 
be parallelized effectively. Calculations reported in this 
work were done on the University of Florida High Per- 
formance Computing Center Linux clusters. 

Though these are quite acceptable costs for fixed ge- 
ometries and small numbers of ions and electrons, the 
burden becomes formidable for direct application in 
Born-Oppenheimer molecular dynamics. We are cur- 
rently working on ways to ameliorate that problem. 
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Appendix A: Correction for Piecewise Integrals for 
the Basis Functions 

From the definition of the basis functions in Eqs. ([7]) 
and ©, it follows that derivatives of the basis function 
may not be continuous at x c , the center of the func- 
tion. A discontinuity of the second derivative would, of 
course, be significant for the kinetic energy. The issue is 
whether the kinetic energy matrix elements can be eval- 
uated piecewise, as is the case with the overlap, nuclear- 
electron, and electron-electron integrals. We may exam- 
ine this issue by writing the basis formally with Heaviside 
functions, as follows: 

Lp = [6{x) - 9{x - xi)] cpo + [0(x - xi) - 6{x - L)] (p L 

(Al) 

For the tGTO basis, the identification from Eq. [S]is 

¥>o = ao [g n {x) - 5 Q ) 
<Pl = a L (g n (x) - S L ) 

g n {x) = (x- Xl ) n e~ a(x ' Xl)2 , xi = x c (A2) 

The derivatives are 
dtp 

— = [S(x) - 5(x - x%)] ip + [0{x) - 9(x - x\)] (p' 

+ [S(x - x\) - S(x - L)] ip L 

+ [6(x-xi)-6(x-L)]tp' L (A3) 

and 

^2 = W( x ) ~ § '( x - x i)] <A) + 2 [S(x) - S(x - xi)] ip' 
+ [8(x) - 9(x - xi)] <p'q 
+ [S'(x - xi) - S'(x - L)] ip L 
+ 2 [5(x — xi) — 5(x — L)] tp' L 
+ [9(x-xi)-9(x- L)]ip'l (A4) 

Now consider a generic kinetic energy matrix element 
involving the foregoing function and another, similar ba- 
sis function \ with left and right constituents xo, XL, 
centered at x^. Without loss of generality, take x\ < X2- 
Then 

f L d 2 w 

J X dx^ dx := Ia + Ib ^ 
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The terms of the second derivative with the Heaviside 
functions contribute just the piecewise integration I a, 
while the delta function and first derivative delta function 
terms contribute Ib ■ Of the terms in Eq. IA4I only those 
at Xi contribute, as the constituents xoi XL, fo, and ifiL 
go to zero at x — and x = L. Thus 

Ib = J Xo [S'(x - xi) {ip L - ipo) 

+28{x - x x ) (ip' L ~ tp' )] dx . (A6) 
From the definition of 5' this expression becomes 



(po) dx 
(A7) 



which reduces to 



Xo [S(x ~ xt) (ip' L - ifo) 
+2S(x - xt) {ip' L - ip' )] - x <5(a; - Xt) (<Pl 



[Xo {<Pl - <Po) - Xo (<Pl - <Po)] 



(A8) 



For the case of xt = x%, that is, for diagonal terms or 
functions that have the same center, x'o must be replaced 
in Eq. ([AT]) by the analog of Eq. (|A3[) . with the result 



Ib = Xo {<Pl - <Po) 



8(x - xt) {<Pl - <Po) [[S(x) - S(x - xt)} xo + [6(x) - 0(x - xt)] x'o 

+ [S(x - xt) - 6(x - L)] xl + [0(x - xt) - 9{x - L)\ X ' L ] dx (A9) 



Note xo m t ne first term follows because the func- 
tions themselves are continuous, while the first and second 
derivatives may not be. For the same reason, it follows 
that the remaining integral in Eq. (IA9|) and the second 
term of Eq. (|A8j) vanish. Thus so long as the functions 
are continuous, the correction to the piecewise kinetic 
energy integral is simply 



X&'l " fo, 



(A10) 



and so long as the first derivative is continuous, this re- 
duces to zero. 

With continuity of the basis functions enforced by con- 
struction, only the first derivative needs to be examined 
for a possible correction to simple piecewise integration. 
For the basis defined in Eqs. (jHJ), those corrections are 



Appendix B: Finite- Range Gaussian Integrals 

Following the methods of Boys^— , we use the trans- 
form of the Coulomb potential to separate the Coulomb 
integrals into one-dimensional Cartesian pieces: 



1 



1 



|r-R N 
1 



e _ s2(r _R N) 2 ds 



Z s 2 (x-X N ) 2 e -s 2 (y-Y N ) 2 e ~s 2 (z~Z N ) 2 dg 



(Bl) 

Hence all Coulomb integrals require integration over the 
transform variable s, which is done by Gauss-Laguerre 
quadrature. 

For ID primitive tGTOs, we note the required finite 
range (a, b) integrals are of the form 



g°(x) = e - a{x - x ^ 

g 1 {x) = {x - Xl )e~ a{x - Xl)2 

g 2 (x) = (x - xtfe-^ x - x ^ 



Ib=0 

Ib = x(o>l ~ do) 
I B =0 (All) 



Basis functions g n (x) with higher powers of the prefactor 
(x — xt) all have continuous first and second derivatives 
at xi. In fact, those derivatives are all zero. So only the 
p-type basis functions, (n = 1), have a kinetic energy ma- 
trix element contribution beyond that given by piecewise 
integration. 



In = I x e 

J a 

This simply transforms to 



n e - a {x-x c f dx 



I n = (x' + x c ) n e- ax dx' 

J a—x c 

This result leaves us needing to compute 
J n = / x n e- ax2 dx 



(B2) 



(B3) 



(B4) 
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Integrating by parts we find 



2a 



n-1 
~2a~' 



>n-2 



(B5) 



So with the initial two integrals, we may find the higher- 
order integrals by recursion: 



dx 



Ji 



-J. 

J a' 



b' 2 ■ " ■ 

xe~ ax dx = - 



e 



erf (yfax) 

b' 



2a 



All two-center (overlap, nuclear-electron, kinetic) inte- 
grals reduce to expressions in terms of /„. After two ap- 
plications of the Gaussian product rule, four-center (two- 
electron) integrals reduce to terms of the form 



b 2 nb 



Two further applications of the product rule bring us to 
the form 



b 2 rb 



c™ x ™ e ~^^- R ^))\- K ^- R *) 2 d Xl dx 2 



(B9) 

(B6) 

Here the integral over x\ may be evaluated as /„, so we 
are left with 

(B7) 



pb 2 

/ I n {X2) 
J a 2 



xi=bi 



m -k 2 (x 2 -R 2 ) 2 



Xn C 



dX2 



(BIO) 



x^x^e- ai{xi - x ' )2 e - a2{x2 - x ^ e - s2 (*i-x2) 2 dxidx2 

(B8) which we evaluate by Gauss-Legendre quadrature. 
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